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Employing simplified models in computer simulation is on the one hand often en- 
forced by computer time limitations but on the other hand it offers insights into 
the molecular properties determining a given physical phenomenon. We employ this 
strategy to the determination of the phase behaviour of quadrupolar fluids, where 
we study the influence of omitting angular degrees of freedom of molecules via an 
effective spherically symmetric potential obtained from a perturbative expansion. 
Comparing the liquid-vapor coexistence curve, vapor pressure at coexistence, inter- 
facial tension between the coexisting phases, etc., as obtained from both the models 
with the full quadrupolar interactions and the (approximate) isotropic interactions, 
we find discrepancies in the critical region to be typically (such as in the case of 
carbon dioxide) of the order of 4%. However, when the Lennard- Jones parameters 
are rescaled such that critical temperatures and critical densities of both models 
coincide with the experimental results, almost perfect agreement between the above- 
mentioned properties of both models is obtained. This result justifies the use of 
isotropic quadrupolar potentials. We present also a detailed comparison of our sim- 
ulations with a combined integral equation/density functional approach and show 
that the latter provides an accurate description except for the vicinity of the critical 
point. 

I. INTRODUCTION 

The study of fluid phase equilibria by computer simulation methods has 
become an extremely active field, since accurate information on thermodynamic properties 
of simple and complex fluids and their mixtures is of enormous importance for a variety 

nun 

of applications |6|, ITJ, [8(. The problem of understanding the phase behavior of such fluids 



is a fundamental problem of statistical mechanics, too 0, [lo|. While such properties in 
principle can be found from experiments, particularly for mixtures such data still are rather 
incomplete, since a cumbersome study of a large space of control parameters (temperature 
T, pressure p, and mole fraction(s) x (x a ) in the case of binary (multicomponent) mixtures) 
needs to be made. 

In simulations an economical use of computational resources often dictates the use of 
models that are as simple as possible. The standard approach for the simulation of fluids 
is to apply classical Monte Carlo and Molecular Dynamics methods 

US 

, that require 

effective potentials (usually of pairwise type). Thus, both quantum effects associated with 
the finite mass of the nuclei are ignored as well as the degrees of freedom of the electrons 
(sometimes the latter are considered, when the effective potentials are derived by "ab initio" 
quantum chemistry methods, see ej 



on purely empirical grounds Q, ll5| ) . 

Now, even when the above approximations are accepted, there the question remains 
whether an all-atom model is needed for the description of intermolecular forces, or whether 
further degrees of freedom may be eliminated. For example, consider carbon dioxide (CO2), 
which is an extremely important fluid due to its use as supercritical solvent 6, Models 



[12I Q , but sometimes these potentials are postulated 



used for tlx 1 simulation for C() 2 arc truly abundant jl^l [la 
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261 ] : they include all-atom models with either flexible or rigid intermolecular distances, 



as well as models where a CO2 molecule is reduced to a point particle, with [27] or even 
without [28[ a quadrupolar moment. While the latter case, where the molecules are described 
as simple point particles interacting with Lennard- Jones forces, is computationally most 
efficient, it is also the least accurate. Recently it was suggested [271, |29( that a significant 



gain in accuracy with almost no loss in computational efficiency can be obtained by using 
perturbation theory to construct an effective isotropic quadrupolar potential 3fj| . While very 
promising results for a variety of molecular fluids, including carbon dioxide and benzene, 27] 
have been obtained, it still needs to be established to what extent the isotropic quadrupolar 
potential actually reproduces the physical effects of the actual angle- dependent quadrupolar 
interactions. 

In the present paper we fill this gap, using the case of CO2 as an archetypical example. 
Applying the same grand-canonical Monte Carlo techniques in conjunction with successive 
umbrella sampling 311] and finite-size scaling analysis 32j, |33( that were used for the work 
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applying the isotropic quadrupolar potential [27|, we obtain for the present model (which is 
described in Sec. 2) the phase diagram in the temperature-density and temperature-pressure 
planes, as well as the interfacial tension (Sec. 3). In Sec. 4 we discuss the optimized choice of 
the Lennard- Jones parameters, and show that the differences between the models with the 
full and averaged quadrupolar interactions are rather minor, when the critical temperatures 
and densities are matched. Sec. 5 presents a comparison with integral equation/density 
functional calculations and shows that the latter approach can describe isotherms as well as 
the coexistence curve of the model very accurately, but is not applicable very near to the 
critical point. Sec. IVII summarizes our conclusions. 



II. MODEL AND SIMULATION METHODS 

A. Models with Full and Averaged Quadrupolar Interaction 

We start from a system of uncharged point particles which have a quadrupole moment 
Q and interact also with Lennard- Jones (LJ) forces, 



U, 



LJ 

ij 
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(1) 



with = \fi — fj\ being the distance between particles i and j at sites r i: fj, and the 
range and strength of the LJ potential are denoted as a and e, respectively. 
The quadrupole-quadrupole interaction is 



with 



34] 



tjQQ = 3_Q_ rQQ 
ij 4 r 5 Jij 

ij 



(2) 



f®. Q = 1-5 cos 2 6,-5 cos 2 6^ + 17 cos 2 6; cos 2 6^- + 2 sin 2 6, sin 2 Oj cos 2 ($, - 

— 16 sin 0j cos 0j sin 9j cos 9j cos($j — $j) , (3) 

where (0,, $j) are the polar angles characterizing the orientation of the uniaxial molecule 
relative to the axis connecting the sites fi,fj of the two particles. 

In order to speed up the Monte Carlo simulation, we wish to introduce a cutoff r c such 
that the total potential is zero for r > r c . This needs to be done such that the total potential 
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is continuous for r = r c , and the same condition should apply for the equivalent isotropic 



quadrupolar potential, resulting from treating the quadrupole-quadrupo 
ond order thermodynamic perturbation theory in the partition function 
we use the following potential ("F" stands for "full potential") 



e interaction in sec- 



29 



30j. As a result, 
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r < r r 



(4) 



, r > r c , 



where Sq = (cr/r c ) 6 — (a/r c ) 12 . For <C r c (and large enough r c ) Eq. (TjJ reduces to 
Eqs. (CE])-©, noting the abbreviation 



Q' 2 



Following previous investigations 



use some results of Ref. 



Qf 

28| we use r c 



(5) 

2- 6 \/2 a. Indeed in this work we will 



271 ]. like the critical lines, which depend on the choice of the cutoff. 



Differences in the thermodynamic properties arising from a different choice of the cut-off are 
minor at least for simple Lennard- Jones potentials (and a suitable renormalization of the 
LJ parameters) [27]. 

Note that the potential in Eqs. ([I])-© is cut off in Eq. (T4J) in such a way that the cutoff 
r c does not depend on the angles 4>i,Qi. The corresponding isotropic spherically averaged, 
potential is ("A" stands for "averaged potential") [29] 



U A 



12 



7 n 



10 

^ +S A 

13 ' 



r < r r 



r > r r 



(6) 



where again the constant Sa is chosen such that the potential is continuous for r 



27] . Note that for the potential, Eq. ([6]), the forces at r c are discontinuous but do not 



diverge there. This is a requirement if one wishes to estimate the pressure from the virial 
theorem, when one does a simulation in the /xVT or NVT ensemble, respectively 

In Ref. J27] we investigate extensively the use of potential for modeling quadrupolar 
substances like carbon dioxide and benzene. These results were compared with prior inves- 



tigations of a simple LJ model 



28 1 without quadrupolar moment in which we only match 
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critical temperature and density with corresponding experimental values to obtain and 
a a- (In principle, and oa could also be matched at any point in the phase diagram, 
which by definition would improve agreement around this point - however, at the cost of 
imposing an inaccurate description of the critical region.) The introduction of a spherically 
averaged quadrupolar moment improves agreement with the experimental phase behaviour 
significantly, especially for carbon dioxide. (Compare with Figs. El [7] and [8]). 

The optimal choice of the parameters sa, ca and qA is somewhat subtle. A straightforward 
choice simply requires that the two potentials are strictly equivalent for temperatures T — > 
oo, where the perturbation expansion becomes exact. This would imply 

e A = e } a A = a } g A = = (7) 

Note that for the potential the parameter is not a constant, but inversely propor- 
tional to temperature T. 

However, the physically most interesting region of the system is clearly not the regime 
T — > oo, but rather the vicinity of the critical temperature T c . Thus, in our previous work 



on carbon dioxide (C0 2 ) [27|] where Eq. fl§D was used, we have chosen the parameters e^, cr A , 
such that both T c and the critical density p c of the model precisely coincide with their 
experimental counterparts [35] T cexp and p Ciexp . Using also the experimental value of the 



quadrupole moment for C0 2 , Q = 4.3 DA, this implies 27] 



e A = 3.491 x i(T 2i J, a A = 3.785 A, q A , c = Qa(T c ) = 0.387, q F = 0.682. (8) 

In Sec. 3, we shall present numerical results for thermodynamic properties of the full 
model, Eq. (j3J) and (j3J) with this choice of parameters, Eq. (jSj), and compare them to the 
corresponding results based upon Eq. (jSJ) (some of the latter results have been compared in 



271 ] to both experiment and simulations of CO2 using other models). 

As we shall see in Sec. 3, in the critical region both models, Eqs. (j4j) (JSJ) and Eq. (jHJ) are 
no longer strictly equivalent to each other, as expected since the accuracy of perturbation 
theory deteriorates the lower the temperature. Being interested in the critical region, it is 
more natural to choose the parameters e, a and €a, cta °f both models such that T c , p c of 
both models match their experimental counterparts. This requires necessarily a choice of e 
and a different from the choice implied by Eqs. (jjTJ), flS}, since the latter choice would yield 
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(choosing Lennard- Jones units e, a of Eq. 

T* {F) = 1.167, p* c {F) = 0.340, (9) 
T* {A) = 1.203, p* c {A) = 0.347, (10) 

as shall be discussed in more detail in Sec. 3. 

Since the relation between q F and Q (Eq. [5]) or q A and Q (Eq. [7J depends on e and 
a as well, finding the choice of the latter parameters which yields T c = T Cjexp and p c = 
Pc, exp is a self consistency problem j^. In principle, one needs to record the functions 



Tc(<1f) /T*(qF = 0) and P*(qf) / Pc(Qf = 0), similarly as described in [27J]. However, since 
the differences between the results of Eqs. fl9]) and ffTO]) are rather small, it is a very good 
approximation to simply keep the value of qp as found in Eq. (jSJ) and just recompute the 
appropriate values of e and a, which we shall denote as ep and a F , in order to distinguish 
them from the choice of Eqs. (|7ll8l . This procedure immediately yields 



ep = 3.598 x 10~ 21 J , (x F = 3.760A . (11) 

A good test of possible errors introduced by this approximation is provided by using 
Eq. ((SD together with Eq. (fTTT) to check the precise value of the physical quadrupole moment 
strength Q this corresponds to. This yields Q new = 4.292 DA instead of the value used in 



27|], Q = 4.300 DA (note that the actual quadrupole moment of CO2 is negative, but since 
only the square of Q actually matters, cf. Eq. (j2J), we suppress the sign of Q throughout). If 
one uses this slightly modified value Q new instead of Q in Eq . ([7 1) together with the master- 



curves T*{q A {T c ))/T*{0) and p*{q A {T c )) / p*(0) calculated in [2j|, instead of Eq. © slightly 
revised estimates of e A and a a would result 

e A = 3.494 x 10~ 21 J , a A = 3.784 A, q A (T c ) = 0.385 . (12) 

But in view of the large error with which the actual quadrupole moment strength of CO2 
is known 35J, Q = 4.3 ± 0.2 DA, Eq. f fl2|) is as good as a representation of reality as the 
choice of [271] (as quoted in Eq. (JE|)) has been. 

In Sec. 4, we shall compare the result of the model with the full quadrupolar interaction 
Eq. (J4j) , using Eq. (fill as LJ parameters, with the averaged interaction, Eq. ([6]), using 
Eq. (|T2|) as choice for the LJ parameters, since then all physical parameters (T c , p c , Q) that 
coincide with their experimental counterparts, have precisely the same values. 
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B. Simulation Methods and Tools for the Analysis of the Simulation Data 



In this section, we summarize our procedures for carrying out the simulations and their 
analysis only very briefly since more detailed descriptions for similar models can be found 



in the literature 



27 



28] 



The estimation of vapor-liquid coexistence curves and critical parameters is done in the 
grand canonical (pVT) ensemble, varying the chemical potential p and recording the density 
distribution Pl{p), and analyzing carefully the dependence on the linear dimension L of the 
simulation box (as usual, we take V = L 3 at the critical point, while deep into the coexistence 
region we use an elongated box V = 2-L 3 . For both geometries periodic boundary conditions 
are applied). For T < T c , the value of p coe x{T) where phase coexistence between vapor 



and liquid occurs is found from the "equal weight rule" |31l . |32| . 1331 . |36| . For an accurate 
sampling of Pl(p) including the densities inside the two phase coexistence region that also 
need to be studied for an accurate estimation of the weights of the vapor and liquid phases, 



31 



33 



successive umbrella sampling methods 23, [31[ are used, as well as re-weighting procedures 



371 ] . Note that the presence of the orientational degrees of freedom in Eqs. ( 131T41) 



does not constitute any principal difficulty here. The acceptance rate for the insertion 
of particles with a randomly chosen orientation is of the same order as for the isotropic 
potential, Eq. ([6]), where this degree of freedom has been eliminated. This fact is understood 
easily, since the strength of the quadrupolar interaction, Eq. ([2]), is distinctly smaller than the 
strength of the Lennard- Jones interaction, Eq. (TjQ), for the present choice of g^. However, the 
time required to compute the energy change caused by such a particle insertion or deletion 
is about an order of magnitude larger when Eq. (BJ rather than Eq. ([6]) is used, due to the 
complicated angular dependence of the quadrupole-quadrupole interaction (Eq. [3]). 

Nevertheless it is still feasible for this model, Eq. (j2J), to obtain sufficiently accurate 
information on Pl{p) for a variety of temperatures T and lattice linear dimensions L, 
following a path along the coexistence curve p = p CO ex(T) in the (p,T) plane, and its 
continuation for T > T c (there the path is defined by the condition that the derivative 
(dp/dp)T = L 3 ((p 2 )T,fj, — (p)r M ) is maximal). Fig. CD shows, as an example, second and 
fourth order cumulants £7 2 and £7 4 along such a path as a function of temperature. These 
cumulants are defined by 
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U 2 = (M')/(\M\y, U 4 = (M 4 )/(My , M — p — (p), 



(13) 



where now we have omitted the subscripts (T,p) from the averages (...). As is well known 



flfl 



11 



32 



331 ]. accurate estimates for T c can be obtained from the common intersection 



point of either U2(L,T) or [7 4 (L,T) for different L. The justification of this simple recipe 



follows from the theory of finite size scaling 



38 



401. Note also that the ordinate values 



11%, Ul of these common crossing points should be universal for all systems belonging to the 
Ising model universality class, to which both models Eq. (TjO) and ([6]) should belong, and are 
known with very good accuracy (3. I^l ] . 

From Fig. [1] it is evident that this intersection property does not work out perfectly 
well, in particular, the curve for L = 6.74 a is somewhat off. However, finite size scaling 
should become exact only in the limit where both L — > 00 and T — > T c , while otherwise 
corrections come into play. Systematic improvements (taking the so-called "field mixing" 
4|] and "pressure mixing" effects [42] into account) are possible, but are not considered 
to be necessary here, since the relative accuracy of our estimate for T c extracted from 
Fig. [I] is clearly not worse than 3-10 -3 , and this suffices amply for our purposes. A recent 
comparative study of different finite size scaling based approaches for the study of critical 
point estimation of Lennard Jones models 43j is in full agreement with this conclusion. Note 
ihat the accuracy of the data in Fig. [T] is comparable to data taken for a pure LJ model 
281 ] or for Eq. (J6j) 27], respectively (we estimate the relative accuracy of the curves in Fig. 
□ to be of the order of 0.5% or better). We have used 7-10 6 , 3-10 6 , 6-10 6 and 9-10 6 Monte 
Carlo steps (respectively for the L/a = 6.74, 9, 11.3 and 13.5 system) for each simulation 
point (T*,A/i coex (T*)), for which the data for -Pl(p) were sampled, and applied histogram 
extrapolation methods (see 



33| for details) to obtain the smooth curves drawn in 
Fig. [H In every step 100 attempts to insert or delete particle plus local moves are done. 

For T < T c the densities pcoL, piolx of the two coexisting vapor and liquid phases can be 
simply read off from the peak positions of Pl{p), and from the density minimum in between 
the peaks the interfacial tension j(T) can be estimated, using the relation 
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28 
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(T)/k B T = 0.5L- 2 HPl(pS£>)/Pl(p*)], 



00 



(14) 
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where 



(pcoL + Pcoex)/2 denotes the density of the "rectilinear diameter". All these 



-,( 2 ) 



3 if ] and their exten- 



methods work well for fluid models with simple isotropic potentials [27|, 
sion to the present model (Eqs. [3j H|) is fairly straightforward. 

For a comparison between the two models defined by Eqs. (j4j) and (EJ) outside of the 
critical region it is also of interest to apply NVT and NpT ensembles. Then no particle 
insertions or deletions occur, but rather a particle is selected at random and a move is 
attempted where one puts it in a randomly chosen position inside a sphere of radius 5r 
around its old position. Simultaneously the orientation of its molecular axis is chosen inside 
a cone of angle 5Q around its old orientation [l, 2|. The choices of 5r and 50 were adjusted 
to have an acceptance rate of 40% for such moves. 



III. DIRECT COMPARISON BETWEEN RESULTS FOR THE FULL AND 
CORRESPONDING AVERAGED QUADRUPOLAR POTENTIAL 

As noted in Sec. 2.1, Eq. (jBJ) results from Eqs. ( 13141) when one carries out a second order 
thermodynamic perturbation theory and interprets the result as being due to an average 



potential 29J. Since thermodynamic perturbation theory is basically a high temperature 
expansion in powers of 1/T, it is a matter of concern how accurate such a procedure really 
is in the temperature region around criticality and below. 

As a first test, we have carried out a NVT simulation using the averaged model, Eq. (JSj), 
at a density that is much larger than the critical density, namely p* = 0.544, and we 
have recorded the corresponding pressure p*(T) from the virial formula (Fig. [2^). This 
pressure then was used as an input for a NpT simulation of the full model, Eqs. (J3J 



Fig. [2](b) shows the corresponding comparison: one sees that for large T* (i.e., T* > 2.5), 
the data obtained from the full potential indeed converge against the density p* that was 
chosen, while for T* < 1.5 there are rather pronounced deviations. Of course, if an NpT 
simulation is carried out using Eq. ([6]), the chosen density p* = 0.544 is reproduced over 
the entire temperature interval shown in Fig. Mb) with negligibly small errors, since with 
the chosen volume [Lja^ = 10.3) systematic discrepancies between the different ensembles 
of statistical mechanics are completely negligible (although such discrepancies will occur in 
the two-phase coexistence region or near the critical point) . The deviations seen in Fig. EJ^b) 
simply represent the higher order terms of the 1/T expansion, by which the averaged and 
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full potentials differ. Similar discrepancies between the full and averaged potential were also 
seen on the vapor side of the coexistence curve (not shown here). 

As a result, coexistence densities in the (T*,p*) plane are rather different as well for the 
two models, Fig. [3j as expected from the differences noted in Eqs. ( l9lfT0l) . and a similar 
discrepancy occurs between the predictions for the pressure along the coexistence curve 
(Fig. HJ) and the interfacial tension (Fig. [5]). 

IV. COMPARISON BETWEEN RESULTS FOR THE FULL AND AVERAGED 
QUADRUPOLAR POTENTIALS WITH OPTIMIZED PARAMETERS 

Now we present a comparison between the full quadrupolar plus Lennard- Jones potential 
(Eqs. HI [5]) and the spherically averaged one (Eq. [6]) choosing the parameters as given in Eq. 
( fTTT) for the full potential and in Eq. ( [12]) for the averaged one, for which critical temperatures 
T c and critical densities p c coincide with their experimental counterparts in both cases. 

Fig. [6] shows that along the vapor branch of the vapor-liquid coexistence curve the av- 
eraged potential slightly underestimates the experimental vapor densities, while the full 
potential slightly overestimate them. However, these deviations are of the same order in 
both cases, and hardly visible (on the scale of Fig. [6]) anyway. Recalling also the fact that 
the coexistence curves (and other data extracted from the simulation) still may suffer from 
systematic effects (residual finite size effects) and statistical errors, see Sec. [Til of the order of 
up to 0.5%, one should not pay too- much attention to the residual differences. We conclude 
that both models describe the vapor branch of the coexistence curve equally well, over the 
studied range of temperatures (which extends from T c down to about T = 250 K). 

However, for the liquid branch of the coexistence curve the model with the averaged 
interaction performs distinctly better. Of course, there is no physical reason known to 
us why this should be the case. We believe that this more accurate description of the 
isotropically averaged model is only due to a fortunate compensation of errors. 

With respect to the vapor pressure at phase coexistence (Fig. [7]), we see, however, that 
at low temperature (250 K < T < 280 K) the model with the full quadrupolar interaction 
performs slightly better than the isotropically averaged model. Near the critical point, 
however, the isotropic quadrupolar interaction performs slightly better, since it predicts 
the critical pressure a bit more accurately. Thus, we conclude that the vapor pressure at 
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coexistence is predicted by both models about equally well. 

Fig. [8] finally compares both model predictions with the data 35| for the surface ten- 
sion between both phases. In this case there is a clear preference for the model with the 
isotropically averaged quadrupolar interaction. Taking the results from Figs. EH7| together, 
we conclude that for a description of phase coexistence the model Eq. ([6]) is clearly the 
better "effective" model. Also for a supercritical isobar (Fig. [10] presents an example) the 
model Eqs. (j4ll5]TTTl) does not have an advantage. The comparisons presented in this section 
thus fully justify the use of Eq. (jSj) for practical purposes. 

An additional interesting test now concerns the temperature dependence of the density 
that results when we compare NVT simulations for the averaged potential with NpT simu- 
lations for the full potential (similarly to what was done in Fig. [2]), choosing the parameters 
of Eq. (|lip and ()12p . We have also included a comparison with two analytical approaches, 
namely an integral equation/density functional theory (IE/DF), see the following section, 



and perturbation theory combined with mean spherical approximation (MSA), see Ref. 27] 



for a description of this method in the current context. Both approaches agree with our 
results very well. Fig. E^a) is the counterpart of Fig. EJ^a); again it is seen that the pressure 
at p = 0.733 g/cm 3 for the full model is in very good agreement with the corresponding 
experimental data and averaged potential. However in this case the full model is superior 
with respect to the averaged model. Fig. Mjo) shows that indeed the NpT results for the 
full potential now converge rapidly to a somewhat higher density (near p « 0.75 g/cm 3 ) 
as the temperature is raised from the critical region to higher temperatures. Of course, 
as expected, it is not possible to fit the critical region (as done in Fig. EHH]) and the high 
temperature region (as done in Figsj2Jl5]) simultaneously. 



V. INTEGRAL EQUATIONS WITH REFERENCE FUNCTION ALS 

In this section we summarize a combined integral equation/density functional method 
to calculate equations of state. A novel approach to avoid unphysical no-solution domains 
near the critical point is outlined and data for the equation of state of the averaged model 
defined by Eq. ([6]) are compared with the simulation results. 

The pair correlation function h(r) and the direct correlation function (of second order) 
c(r), which contain all thermodynamic information of a given homogeneous model system of 
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density P n an d interacting with an isotropic pair potential £/y(r), fulfill the following general 



relations 



44 



45 



46|: 



h{r) — c{r) = po j d 3 r' h(\r - r'|) c(r') , 
ln[l + h(r)\ - pUij(r) = h(r) - c(r) - b(r) . 



(15) 
(16) 



The first relation is known as the Ornstein-Zernike equation, while the second is the general 
closure relation in terms of a yet unknown bridge function b(r). For a solution, it is necessary 
to specify the bridge function b in terms of h and c. Thermodynamics is obtained either 
through the virial route [47J, giving the pressure p, 



^ = 1- np (3 rfrr 3 (l + /i(r))— ^ 
Po 3 Jo dr 



or through the compressibility route 47] 

dp 







dpo 



An po I drr 2 c(r) . 
'o 



(17) 



(18) 



Both routes are not necessarily identical for a given approximation for the bridge function. 

For one-component systems, advanced methods exist which yield good agreement with 
simulations for the pressure in the whole p—T plane and for the whole coexistence curve, e.g. 
the SCOZA (self-consistent Ornstein-Zernike approximation) 48] or the HRT (hierarchical 
reference theory) 49] . A drawback of these methods is their "non-locality" in the p — T 
plane, i.e. for a SCOZA solution a partial differential equation has to be solved on this 
plane and for a HRT solution a renormalization flow equation has to be solved on a specified 
isotherm. 

Therefore we treat our problem within a formalism which is close to the reference 
hypernetted-chain (RHNC) method. In its original version 50] (developed for repulsive 
core fluids), the bridge function b was taken from a reference hard sphere system with 
suitably optimized reference packing fraction (or hard sphere diameter). Since the RHNC 
closure equation ([TBI) can be derived from an approximate bulk free energy functional, a 
closed expression for the chemical potential /i also exists. Near the critical point, however, 
there exists a no-solution domain (extending into the supercritical region T > T c ) where 
no real solution to the RHNC closure can be found. As shown below, this problem can be 
eliminated by a method (FHNC) where a generating functional for the bridge function is 
adopted from a suitable reference system. 
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A. Bridge functions from a reference functional 



The FHNC generalization of RHNC within the framework of density functional theory 
was proposed in Ref. 5l| . By making a Taylor expansion of the free energy functional 
around bulk densities the general closure equation is derived and the bridge function in the 
bulk is determined via a density functional for a suitably chosen reference system of hard 
spheres. To be more explicit, the excess free energy functional 



^ x [p(r)]=^ HNC [p(r)]+^>(r)] 



(19) 



is split into a part (^ rIiNC ) which generates the hypernetted-chain (HNC) closure {b = 0) 
and a remainder (J 73 ) which generates a non-zero bridge function. The HNC part jF HNC is 
a Taylor expansion of the excess free energy up to second order in the deviations from the 
bulk density po, Ap(r) = p(r) — po: 

jtHnc = F ( po) + ^cx J rf s r Ap(r) _ J_ J d s r j rf 3 r / c( | r _ r'|) Ap(r) Ap(r') . (20) 

There the defining relations for the excess chemical potential p ex (po) and the direct correla- 
tion function c have been used: 



5F C 



5p(r) 



A*"(pd) 



2 tcx 







6 2 F 



-c(|r - r' |; po) • 



(21) 



p(r)=p(r')=po 



p(r)=p () 5p(r)5p(r>) 
The general closure equation follows by employing the test particle method: the grand 
potential 



Q[p(r)} = ^ d [p(r)]+^ x [p(r)]- / rf»r (p - V(r))p(r) 



f3r d [p(r)] = J d 3 rp(r)(ln[A 3 p(r)] - 1) 



(22) 



is minimized in the presence of a fixed test particle of the same species which exerts the 
external potential V(r) = Uij(r) onto the fluid particles 5l|, [52| (A is the thermal de-Broglie 
wavelength). The precise form of the closure equation (Eq. [16]), is recovered upon the 
following identifications: 



h{r) 



Ap(r) 



6(r) = /3 



5F l 



(23) 



p(r)=p (h(r)+l) 



Po Sp(r) 

In general, the excess free energy functional beyond second order, the bridge functional J-" B , 
is not known. Therefore the key step of the present method is to approximate JF B by a 
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density functional for a reference system in the following manner: 

f B [p] « F B > rci [p] = r ci [p] - F UNC ' rei [p] , (24) 



where the second order HNC contribution is defined as in Eq. (12011 with the replacements 
F(po) — > -F rcf (po), p cx —> /i ex,ref and c — > c ref . For fluids with repulsive cores, the reference 
functionals of choice are hard-sphere functionals which are known with high accuracy (such 
as e.g. the ones in Refs. [53, 54]). In such a manner the system of equations (fT5|) and (1TB]) is 

n, 

closed and amenable to numerical treatment. According to Ref. [52| the optimal reference 
hard sphere packing fraction ?7 ref = (ir/6) Pod^ ef is determined through the (local) condition 

-|- [p {h{r) + 1); d rc{ ] - F B ^[p (h^(r) + 1); d Te{ }) = , (25) 

which corresponds to extremizing the free energy difference between the fluid and the refer- 

potential 



52]. Thus 



ence system with respect to the reference hard sphere diameter d Te f. The chemica 
of the fluid can also be expressed as a functional of h locally in the p — T plane 
a coexistence curve for a given fluid can be determined straightforwardly by the equality of 
p and p on the fluid and the gas side, respectively, and no thermodynamic integrations are 
necessary. 



B. The critical region 

Similarly to RHNC and HNC, the FHNC method outlined above, together with the 

optimization criterion for r/ ref , Eq. (|25p . exhibits a no-solution domain in the p — T plane 

which stretches into the supercritical region (T > T c ). This can be attributed to a failure of 

the optimization criterion in the critical region which assigns a wrong long-range behavior to 

the direct correlation function c. Consider the asymptotic expansion of the closure, Eq. (fT6j) . 

where h, c and b are small: 

h 2 (r) h 3 (r) , . T/ , . , , . 

y- + —-^ + ■■■ = -c(r) - b(r) (r -> 00) . (26) 

2 3 

(Here we assume that the potential C/y(r) is cut off.) In HNC (RHNC) the bridge function 
b is zero (short-ranged), therefore we find to leading order c(r) ~ h 2 (r)/2. However this is 
inconsistent with the critical behavior of the correlation functions. In three dimensions, this 



critical behaviour is approximately described by the Ornstein-Zernike form 55 1 
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Through the Ornstein-Zernike equation (fT5l) it follows that in this limit c oz (r) stays short 
ranged in the sense that its Fourier transform c (q) = Cq + c 2 q 2 + • • • permits an expansion 
around q = 0. In Eq. ( |27j) . £ is the correlation length which goes to infinity upon approaching 
the critical point and is related to c oz through £ 2 = — poc 2 /(l — poQ))- Assuming h — > h oz , 
the asymptotic HNC (RHNC) closure c ~ h 2 /2 is in conflict with the requirement of c 
staying short-ranged upon approaching the critical point, i.e. it cannot be a solution of the 
Ornstein-Zernike equation. On the other hand, within FHNC the bridge function b itself 
depends on h as 

r) Ci2 gx rcf 

b(r) » W^Y^y h\r) + 0(h*) (r - oo) . (28) 
and the asymptotic closure, Eq. (|26|) reads 

i X *~\2 gx rcf \ 

< r ) = 2 y 1 ~ ^ )2 ^d(^) ^ + ° {Ki) ■ (29) 

Thus we see that upon requiring /3(r] ref ) 2 g^ rC f)2 = 1 the closure is consistent with h — > h oz 
and c staying short-ranged. This condition is fulfilled for rj mf w 0.13 and in the critical 
region, this condition on the reference system packing fraction replaces Eq. (1251) . Incidentally, 
this condition is consistent with the intuition that the reference hard sphere diameter <i re f is 
roughly equal to the Lennard- Jones diameter a for densities close to the critical density. We 
checked for various supercritical isotherms that the modified optimization criterion indeed 
removes the no-solution domains. (For a different approach to this problem within FHNC, 
see Ref. 13 •) 

C. Numerical results 

Numerical data for the coexistence curve (virial route) of the averaged model (Eq. [6]) 
are given in Fig. [3] (g^ = 0.387) and Fig. [6] (g^ = 0.385). The overall agreement with the 
simulation results is good, except for temperatures within 5% of the critical temperature 
T c where a noticeable shift of the coexisting gas densities towards higher values can be 
observed. This is also the reason why the pressure on the coexistence curve is somewhat 
larger than the pressure determined by the simulations (Fig. [7]). Analyzing the behavior of 
the solutions in the critical region more closely, we find that there (as many other integral 
equation approaches) the FHNC method suffers from the inconsistency between the virial 
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and the compressibility route to the equation of state. Indeed, the compressibility route (for 
qA = 0.387) gives a critical temperature of T c com « 1.14 and critical density p* « 0.36 with a 
behavior of the correlation function h consistent with the Ornstein-Zernike form, Eq. (1271) . 
confirming the applicability of the reasoning in the previous subsection. Therefore, the virial 
route coexistence data for temperatures above T c com are not reliable. 

We remark that the FHNC method is not particularly designed to capture the correct 
critical behavior. With respect to the prediction of the coexistence curve and coexistence 
pressure only, it is not particularly superior in accuracy to the (simpler and faster) per- 
turbative Mean Spherical Approximation (MSA) whose results have been discussed in j^l 
(for a detail description of the Equation of State used in j^, we refer to 57, 58]). Clearly, 
in this respect renormalization-group based methods such as HRT perform much better. 
However, supercritical properties of CO2 are reproduced quite accurately as the comparison 
of the p — p isotherms (T = 316.36 K) between the experimental data for CO2, FHNC and 
MSA shows (Fig. [TTj) . While the experimental data and the FHNC results are almost on 
top of each other, the perturbative MSA results exhibit a van-der-Waals loop due to the 
underlying mean-field approximation which results in a too large T c j^]. Additionally we 
observe very good agreement with simulations for the supercritical isochore p = 0.733 g/cm 3 
(Fig. ED and the supercritical isobars p=100 bar and 200 bar (Fig. [TO]) . 

A problem of the FHNC approach seems to be the accurate prediction of surface tensions. 
Although the technique can be extended to compute this quantity, the results are much less 
satisfactory, since the simulation results are about 40% lower than the FHNC results, in 
the temperature region around T=270 K where the coexistence curve is predicted rather 
satisfactorily by FHNC (Fig. [6]). A similar problem was noted in a recent comparison of 
Monte Carlo and density functional theory results for phase separation in colloid-polymer 
mixtures 59 1. 



In concluding this section, we find that the presented FHNC method allows a fast and 
precise determination of the equation of state except for the vicinity of the critical point. 
Within FHNC, the pressure p and the chemical potential p are obtained through local re- 
lations in the T — p plane. It appears as an advantage that FHNC is straightforwardly 
generalizable to mixtures since the functionals for the reference hard-sphere mixtures are 
known. First studies of Lennard- Jones mixtures [so] confirm the accuracy of the approach. 
Besides the computation of the pair structure in fluids, the connection to density func- 
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tional theory makes FHNC also a versatile tool to study wetting /drying phenomena 52] 



and effective depletion potentials in dilute colloidal solutions 



61 



62 



63. 



64). 



VI. CONCLUSION 

In the present work, two coarse-grained models of quadrupolar fluids were compared 
with each other (and with both experiment and a theory (IE/DF) combining an integral 
equation approach with density functional theory). The aim of our works is to develop some 
understanding for which region of parameters such coarse grained models are accurate, and 
also to clarify the applicability of the analytic theory. While we use experimental data for 
carbon dioxide as a prototype example of a quadrupolar fluid for comparison, our aim is not 
to provide a chemically realistic modeling description of this substance or any other material. 



Recalling that there is a 
for chemical processing [6 



ot of interest to use supercritical carbon dioxide as a solvent and 
SI], and that other quadrupolar fluids such as benzene also find 
widespread applications, there is a need for efficient coarse-grained models of such simple 
molecular fluids. (A chemically realistic modeling of systems like polystyrene-carbon dioxide 
mixtures would be far beyond reach). Due to the fact that critical fluctuations invalidate 
simple analytic theories extending the Van-der- Waals approach (see 65J for a discussion), 
a simulation approach as presented here is well-suited to include a sufficiently accurate 
description of the critical region. 

We model the quadrupolar fluid by spherical point particles carrying a quadrupole mo- 
ment (the strength of this quadrupole moment being taken from experiment), so that the 
total interaction between two molecules is the sum of a Lennard- Jones interaction (Eq. |T|) 
and the quadrupole-quadrupole interaction (Eqs. [2j [3]). Possible three-body forces are not 
at all included explicitly, but to some extent implicitly, since the Lennard- Jones parameters 
of our effective potential are chosen such that the actual critical temperature and density of 
the material (CO2 in the chosen example) are reproduced. For the sake of computational 
efficency, the potential is truncated at a cutoff distance r c and shifted to zero there (Eq. HJ). 

As a second model we chose a closely related one, where the angular dependence of the 
quadrupolar interaction is averaged perturbatively, Eq. (|6|). This isotropic potential can as 
easily be treated numerically (sec. |V|) as other simple isotropic pairwise potentials. 

By construction, the two models have to agree at very high temperatures, but this is not 
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the region of interest for applications. In the critical region, discrepancies of the order of 
4% are typically found, in the case of carbon dioxide. 

However, when we determine the effective Lennard- Jones parameters for both models 
such that they reproduce the same critical temperature and density (namely the critical 
parameters of carbon dioxide in our case), we find that both models give a similarly accurate 
description of the equation of state over a rather wide region of parameters (Fig. [6j [7J 
ITUl) . Considering also the surface tension (Fig. [8]), the simpler model with the averaged 
interactions is in better agreement with experiments, despite the fact that at high enough 
temperatures, relative deviations of the averaged model from experiment of the order of 
1-2% can be identified (see Fig. [9] for example). But even for quantities such as isobars at 
p=100 bar and 200 bar (Fig. [TOi) . about three times the critical pressure, the experimental 
results are very well reproduced over the full density region of interest (from 0.1 g/cm 3 
to 1.0 g/cm 3 ). For such data outside of the critical region, our IE/DF theory yields a 
quantitatively accurate description without any adjustable parameter whatsoever, provided 
we use the Lennard- Jones parameters obtained from the Monte Carlo study as an input 
(fitting the Lennard- Jones parameter to the critical parameters of the analytic theory is not 
appropriate, of course, since the latter is inaccurate). 

Of course the fact that the simple isotropic model is even slightly "better" than the 
more complicated one, as far as the comparison with experiment goes, must be attributed 
to some lucky cancellation of errors. In particular, the temperature dependence of the 
interface tension (Fig. [8]) suggests that the full model might itself be too simple for an 
optimal description of the fluid, and it may be necessary to include additional effects like 
the non-spherical shape of the molecule. (The model with the angular-dependent interaction 
is still far from a full description of chemical reality, of course: but the comparison presented 
in 27]] shows that many very sophisticated atomistic models do not perform better than the 
current simple isotropic model either). A possible explanation of this could stay in the fact 
that the physical quadrupolar moment used in coarse grained interactions (Eqs. H] and [6]) 
could require some effective corrections related to the truncation of our potentials. However, 
also in view of our need of efficient coarse grained models, the fact that the averaged model 
is definitively faster than the model in which angular degrees of freedom are taken into 



account, leads to the cone 



of the averaged models (27J, |29 |. 



usions of the present work, namely we strongly support the use 
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We hence suggest that the present approach using effective isotropic models for 
quadrupolar fluids in spite of the angular dependence of the interactions in such systems is 
useful and we plan to extend it to binary fluid mixtures as well. 
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FIG. 1: Second and fourth order cumulants plotted vs. T* = T/e for qp = 0.682 using 
the model defined by Eqs. fl3|4|) . for four choices of L, namely L/a = 6.74,9, 11.3 and 13.5 
(note that the slope of these curves increases with L). The dotted horizontal lines indicate 
the theoretical values for the 3d-Ising universality class J4I]. 



FIG. 2: (a) Reduced pressure p*(T) (in units of the parameters of Eq. plotted vs. 
reduced temperature T* for the reduced density p* = 0.544, as obtained from the averaged 
potential, Eq. (jSj) (see x). (b) Reduced density p* plotted versus reduced temperature, 
when one takes the pressure from part (a), as input for an NpT simulation, using the full 
potential, Eqs. (jHJ), (BJ, with the parameters chosen as given in Eq. (jHJ) (see o). Note that 
the statistical errors are estimated not to exceed the size of the symbols. 



FIG. 3: Vapor-liquid coexistence curve in the (T*,p*) plane as predicted by Eq. (J6]) 
(full line), using the parameters as quoted in Eq. (jHJ), and as predicted by Eqs. (jHJ), (EJ 
(broken line), using the corresponding parameters (Eqs. E]) implying exact agreement 
between both models in the limit T — > 00. For the averaged model, we also report results 
of the integral equation/density functional theory described in Sec. |V] (see *). The relative 
accuracy of the curves representing simulation results in this figure and in the following is 
estimated to be better than 0.5%. 



FIG. 4: Vapor-liquid coexistence curve in the (p*,T*) plane, for the same choices as 
in Fig. [31 For the averaged model, we report also results of the integral equation/density 
functional theory described in Sec. [V] (see *). 

FIG. 5: Interfacial tension plotted vs. T*, for the two models as specified in Fig. [3j 



FIG. 6: Same as Fig. [3], but choosing the parameters of Eq. (1111) for the model with 
full quadrupolar interaction (broken line) and of Eq. ffl2l for the model with the averaged 
interaction (full line). Experimental data from Ref. 35j are included (broken-dotted line). 
With respect to Fig. [3j critical temperature and density for both models now coincide 
with the experimental values. We also include the LJ predictions of Ref. 28j] (dotted line). 
We notice that the spherical and averaged model in this reparametrized plot produces 
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coexistence densities that are in good agreement. Indeed differences are comparable to 



the two models used to describe C0 2 in our previous paper 



27| (i.e. g A , c = 0.387 and 



Qa,c — 0.470) and with discrepancies of the Lennard Jones model in predicting the phase 



diagram for noble gases (see Fig. 11 of Ref. 23]). Finally we note a better agreement of the 



averaged model with experimental results. For the averaged model, we also report results 
of the integral equation/ density functional theory described in Sec. [V] (see *). 

FIG. 7: Coexistence pressure. We report the results for the full model with simula- 
tion parameters given in Eq. ( fill (broken line), the results for the averaged model with 



simulation parameters reported m Eq. « (full Hue) aud the experimental results M 
(broken-dotted line). We also include the LJ predictions of Ref. [28] (dotted line). We 
observe that at low temperatures the full model gives better results with respect to the 
averaged model. This is due to the fact that for high densities the orientational part of 
the quadrupolar interaction becomes more important. On the other hand, near the critical 
point the averaged model performs better than the full model. For the averaged model, we 
report also results of the integral equation/density functional theory described in Sec. fVl 
(see *). 



FIG. 8: Prediction of interface tension. We report the results for the full model 
with simulation parameters given in Eq. (TTT]) (broken line), the results for the averaged 
model with simulation parameters reported in Eq. ( [121) (full line) and the experimental 



results [35[] (broken-dotted line). We also include the LJ predictions of Ref. [281] (dotted 
line). The averaged model is in perfect agreement with experimental results. 

FIG. 9: (a) Supercritical isochore (p=0.733 g/cm 3 ) for the averaged model with pa- 
rameters given in Eq. ( TTT!) (full line) and for the full model with parameters given in Eq. 
(112p (broken line). Experimental results are also included |35] (broken-dotted line). We 
observe a very nice agreement between the experimental results and the full potential. The 
small discrepancy between the averaged model and the full model at high temperatures 
can be understood in the light of the results of Fig. [21 Indeed, due to the fact that the 
use of the same e and a for the averaged and full models (see Eq. [8]) produces the same 
equilibrium states at high temperature (see Fig. EJ), the new choice of parameters Eqs. ([Til 
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( TT2|) produces a systematic discrepancy at high temperature. For the averaged model, we 
also report results of the integral equation/density functional theory described in Sec. [V] (*) 
and the results of the perturbative MSA theory described in [^j (gray line). The inserted 
picture shows as both the theory are in very nice agreement with the MC results, 
(b). Similarly as in Fig. EJ^b) we report the prediction for the densities for the full model 
(diamond) obtained in an NpT simulation taking as input the pressures plotted in Fig. [9[a) 
for the averaged model (full line). The resulting densities agree well with experiments [35]. 
The small deviation at high temperature (between the two models) can be explained by 
considerations discussed for Figs. [2]and[9](a). 



FIG. 10: Supercritical isobars (p=100 and p=200 bar). We report the results for 
the full model with simulation parameters given in Eq. (fTTj) . the results for the averaged 
model with simulation parameters reported in Eq. (Tl2|) and the experimental results 35]. 
The agreement is very good. For a comparison with two other averaged models we refer 
to Fig. 8 of our previous paper j^. For the averaged model, we report also results of the 
integral equation/density functional theory described in Sec. [V] 



FIG. 11: Comparison between the Integral Equation/Density Functional (IE/DFT) 
theory (*) and an e quat ion of state in the perturbative Mean Spherical Approximation 



(MSA) described in 27|, |57|, |58[ (full line) for the supercritical isotherm T=316.36 K. 



IE/DFT, if compared to experiments 35J|, performs better for intermediate densities 0.3 
g/cm 3 < p <0.8 g/cm 3 . However for p > 0.8 g/cm 3 the two theory predict almost the same 
equilibria states. 
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